use rand::Rng;

const MILLER_TEST_COUNT : usize = 12;

fn add_with_mod(a: u128, b: u128, modulus: u128) -> u128 {
    // 带模加法
    let c = a.wrapping_add(b);
    if c < modulus {
        c
    } else {
        c.wrapping_sub(modulus)
    }
}

fn mul_with_mod(mut a: u128,mut b: u128, modulus: u128) -> u128 {
    // 带模乘法
    let mut res = 0;
    a %= modulus;
    while b > 0 {
        if b % 2 == 1 {
            res = add_with_mod(res, a, modulus);
        }
        a = (a << 1) % modulus;
        b >>= 1;
    }
    res
}


fn exp_with_mod(mut base: u128, mut exp: u128, modulus: u128) -> u128 {
    // 带模求幂
    let mut result = 1;
    base %= modulus;
    while exp > 0 {
        if exp % 2 == 1 {
            result = mul_with_mod(result, base, modulus);
        }
        base = mul_with_mod(base, base, modulus);
        exp >>= 1;
    }
    result
}

pub fn miller_rabin_test(x: u128) -> bool {
    // 素数判定
    if x < 2 {
        return false;
    }
    if x <= 3 {
        return true;
    }
    if x % 2 == 0 {
        return false;
    }
    let mut d = x - 1;
    let mut r: u32 = 0;
    while d % 2 == 0 {
        r += 1;
        d >>= 1;
    }
    let mut rng = rand::thread_rng();
    for _ in 0..MILLER_TEST_COUNT {
        let a = rng.gen_range(2..=x - 1);
        let mut t = exp_with_mod(a, d, x);
        if t == 1 || t == x - 1 {
            continue;
        }
        let mut composite = true;
        for _ in 0..(r - 1) {
            t = mul_with_mod(t, t, x);
            if t == x - 1 {
                composite = false;
                break;
            }
        }
        if composite {
            return false;
        }
    }
    true
}

fn gcd(a: u128, b: u128) -> u128 {
    // 最大公因数
    if b == 0 {
        a 
    } else if b == 1 || a == 1 {
        1
    } else {
        gcd(b, a % b)
    }
}

fn pollard_rho(x: u128) -> u128 {
    // Pollard-Rho算法分解因数
    let mut s = 0;
    let mut t = 0;
    let mut rng = rand::thread_rng();
    let c = rng.gen_range(1..=(x - 1));
    let mut val = 1;
    let mut goal: u128 = 1;
    loop {
        for step in 1..=goal {
            t = add_with_mod(mul_with_mod(t, t, x), c, x);
            val = mul_with_mod(val, if t > s { t - s } else { s - t }, x);
            if step % 127 == 0 {
                let d = gcd(val, x);
                if d > 1 {
                    return d;
                }
            }
        }
        let d = gcd(val, x);
        if d > 1 {
            return d;
        }
        goal <<= 1;
        s = t;
        val = 1;
    }
}

fn factor(mut x: u128, max_factor: &mut u128) {
    if x <= *max_factor {
        return;
    }
    if miller_rabin_test(x) {
        *max_factor = x.max(*max_factor);
        return;
    }
    let mut p = x;
    while p >= x {
        p = pollard_rho(x);
    }
    while (x % p) == 0 {
        x /= p;
    }
    factor(x, max_factor);
    factor(p, max_factor);
} 

pub fn find_max_prime_factor(number: u128) -> u128 {
    let mut max_prime_factor: u128 = 1;
    factor(number, &mut max_prime_factor);
    max_prime_factor
}
